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ABSTRACT 

Alfven waves are a particular class of magnetohydrodynamic waves relevant 
in many astrophysical and laboratory plasmas. In partially ionized plasmas the 
dynamics of Alfven waves is affected by the interaction between ionized and neu- 
tral species. Here we study Alfven waves in a partially ionized plasma from the 
theoretical point of view using the two-fluid description. We consider that the 
plasma is composed of an ion-electron fluid and a neutral fluid, which interact by 
means of particle collisions. To keep our investigation as general as possible we 
take the neutral-ion collision frequency and the ionization degree as free parame- 
ters. First, we perform a normal mode analysis. We find the modification due to 
neutral-ion collisions of the wave frequencies and study the temporal and spatial 
attenuation of the waves. In addition, we discuss the presence of cut-off values 
of the wavelength that constrain the existence of oscillatory standing waves in 
weakly ionized plasmas. Later, we go beyond the normal mode approach and 
solve the initial-value problem in order to study the time- dependent evolution 
of the wave perturbations in the two fluids. An application to Alfven waves in 
the low solar atmospheric plasma is performed and the implication of partial 
ionization for the energy flux is discussed. 

Subject headings: Magnetic fields - Magnetohydrodynamics (MHD) - Plasmas - 
Sun: atmosphere - Sun: oscillations - Waves 

1. Introduction 



Alfven waves are a particu lar class of magnetohydrodynamic (MHD) waves driven by 
magnetic tension ( Alfven! 1942 ). In a uniform and infinite plasma the motions of Alfven 
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waves are inco mpressible and polarized perpendicularly to the directi on of the magnetic 



field 



see, e.g., 



Hasegawa fc Uberoil Il982t I Cramer! l200ll ; iGoossend 120031). Alfven waves a re 



found in both laboratory and astrophysical plasmas (see review by lGekelman et al.ll201lf ) 



Since the pioneering works by, e.g., iPiddingtonl (119561 ) and Kulsrud fc Pearcd (119691 ) it 
is known that partial ionization of the plasma affects the dynamics of Alfven waves. The 
feature most extensively investigated in the literature is the wave damping due to collisions 
between ions and neutrals, although other effects as, e.g., the existence of cut-off values of 
the w avelength are also an important consequence of partial ionization (IKulsrud fc Pearce 
1969h . 



Most of the works that studied Alfven waves in partially ionized plasmas adopted the so- 
called single-fluid approximation (see, e.g.. lBraginskiilll965l ). The single-fluid approximation 
assumes a strong coupling between ions and neutrals. For an MHD wave, this condition 
means that the wave frequency has to be much lower than the frequency at which ion and 
neutrals collide. In other words, it is necessary that in one period there are enough collisions 
for ions and neutrals to behave as one fluid. This restriction is fulfilled in, e.g., the partially 
ionized so lar plasma and so the s ingle-fluid approximation is usually adopted i n that case 



see, e.g 



De Pontieu et al. 2001; Khodachenko et al. 2004; Forteza et al. 



2007; Soler et al. 



20091 . among others). 



An alternative approach is the multi- fluid theory (see, e.g., IZaqarashvili et al.ll2011bl ). 
which considers the various species in the plasma as separate fluids. In the multi-fluid 
description no restriction is imposed on the relative values of the wave frequency and the 
collision frequency, although the mathematical treatment is usually more complicated than in 
the single-fluid case. However the multi-fluid theory has the advantage that it is more general 
than the single-fluid approximation and can be used regardless the value of the collision 
frequency. Hence the multi-fluid theory is adequate to study MHD waves in those situations 
where the single- fluid approximation does not apply. This may be th e case of molecular 



clouds (see, e.g., |Pudritzlll990t lBalsaralll996l ; iMouschovias et al.ll201ll ). A particular form 



of the multi-fluid theory is the two-fluid theory in which ions and electrons are considered 
together as an ion-electron fluid, while neutrals form another fluid that interacts with the ion- 
electron fluid by means of collisions. For th e investigation of MHD wav e s this approach was 



follow e d by, e.g., iKumar fc Roberts! ( 120031 ); IZaqarashvili et al.l (j2011bl ); IMouschovias et al. 
d201lh : ISoler et al.l fcp!2h ~ 



Despite the existing literature on this topic (see the references in the above paragraphs), 
the purpose of the present article is to revisit the theoretical investigation of Alfven waves in 
partially ionized plasmas using the two-fluid theory. Our reasons for tackling this task are 
the following. 
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First of all, the existing papers i n the literature often foc us on very specific situations 
as, e.g., the interstellar medium 



Pudritz 



2003 



1990 



Balsara 



1996 



Zaqarashvili et al. 



(e.g., 



Mouschovias et al 



2011b: Soler et al 



Ku lsrud fc Pearcdll969l ). molecul ar clouds (see, e.g., 



201 ll ) , and solar plasmas ( IKumar &: Roberts 



20121 ) among other cases. Here our aim is to keep 



the investigation as general as possible. To do so we take the neutral-ion collision frequency 
and the plasma ionization degree as free parameters. This makes the results of the present 
article to be widely applicable. We put emphasis on the mathematical transparency and on 
the finding of approximate analytic solutions. 

In addition, there has been recently some confusion ab out the existence of cut-off wave - 
length s for Alfven waves i n a p artially ionized plasma (see IZaqarashvili et al.ll2011bl . 120121 ) . 
While IZaqarashvili et al.l ( 120121 ) have shown that the presence of cut-offs in the single-fluid 
approximation is a mathematical artifact, the existence of cut-offs in the two-fluid case is 
a rea l physical phenomenon (e .g., iKulsrud fc Pearcd 1 19691 ; iPudritzl Il990l ; iKamaya fc Nishi 



19981 ; iMouschovias et al.l 1201 ll ) . An important goal of the present article is to stress the 



existence of physical cut-offs for Alfven waves in a two-fluid plasma. 

Finally, unlike previous works that are restricted to the normal mode analysis, here 
we combine results of normal modes with the solution of the initial-value problem. This 
procedure allows us to investigate how strong is the coupling between the perturbations in 
the ionized fluid and the neutral fluid depending on the relative values of the wave frequency 
and the neutral-ion collision frequency. 

This paper is organized as follows. Section [2] contains the description of the equilibrium 
configuration and the basic equations of the two-fluid theory. Alfven waves are investigated 
following a normal mode analysis in Section [3j while the initial- value problem is solved in 
Section HI Later, Section contains an application to the solar atmospheric plasma and 
Section [6] discusses the implications of partial ionization for the energy flux of Alfven waves. 
Finally the conclusions of this work are given in Section [71 



2. Equilibrium and basic equations 



We consider a partially ionized medium composed of ions, electrons, and neutrals. We 
use the two-fluid theory in which ions and electrons are considered together as an ion-electron 
fluid, i.e., the ionized fluid, while neu trals form another fluid that interacts with the ionized 
fluid by means of collisions (see, e.g., IZaqarashvili et al.l 1201 lbl ; ISoler et al.ll2012l ). In all the 
following expressions, the subscripts T and V refer to the ionized fluid and the neutral fluid, 
respectively. 
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The equilibrium is made of a uniform and unbounded partially ionized plasma. We use 
Cartesian coordinates. The equilibrium magnetic field is straight and constant along the 
z-direction, namely B = B z. We ignore the effect of gravity. We also assume that the 
equilibrium is static so that there are no equilibrium flows. The governing equations for the 
various species composing the plasma can be found in, e.g., IZaqarashvili et al.l (|2011bJ). Here 
we restrict ourselves to the study of linear perturbations superimposed on the equilibrium 
state. Hence, the governing equations are linearized. The resulting equations are 



Pi 



Pn 



<9V; 

dt 

dt 
db 

~dt 
dpi 

at 

dt 



-Vpi + - (V x b) x B - on 



-Vpn - « in (V r 
V X ( Vi X B) , 

- 7 PiV • Vi , 
-7P n V • v n , 



Vi 



(2) 
(3) 
(4) 
(5) 



where v i; p- u Pi and p x are the velocity perturbation, pressure perturbation, equilibrium 
pressure, and equilibrium density of the ionized fluid, v n , p n , P n , and p n are the respective 
quantities but for the neutral fluid, b is the magnetic field perturbation, /i is the magnetic 
permeability, 7 is the adiabatic index, and a- m is the ion- neutral friction coefficient. In the 



specific case of a hydrogen plasma, the expression of «i n is given by lBraginskiil (119651 ). namely 



lfhPn / l6k B T 
2 m n V 



(6) 



where m\ and m n are the ion and neutral masses, respectively (mi w m n for hydrogen), k-g 
is Boltzmann's constant, T is the plasma temperature, and <7i n is the collision cross section. 
In the following analysis we do not use this expression of since we take free 
parameter. We do so to conveniently control the strength of the ion-neutral friction force. 
Equation ([6]) is used in the application to solar plasmas done in Section El 

We perform a Fourier analysis of the perturbations in space. In linear theory, an arbi- 
trary perturbation can be represented by the superposition of Fourier components, so that 
we can restrict ourselves to study particular Fourier components. Therefore, the spatial de- 
pendence of perturbations is put proportional to exp (ik x x + ik y y + ik z z), where k x , k y , and 
k z are the components of the wavenumber in the x-, y-, and z-directions, respectively. In a 
uniform and infinite pla sma Alfyen waves are the only MHD modes that propagate vorticity 
perturbations (see, e.g.. ICramerll200ll ; lGoossensll2003l ). In addition, Alfven waves are incom- 
pressible and their motions are confined to perpendicular planes to the magnetic field, i.e., 
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Vi, z = v n,z = 0. Therefore, an appropriate quantity to describe Alfven waves is the vorticity 
component along the magnetic field direction. By working with vorticity perturbations we 
are able to decouple Alfven waves from magnetoacoustic waves. We define 1^ and T n as the 
^-components of vorticity of the ionized fluid and the neutral fluid, respectively, 

Ti = (V x Vi) • z = ik x v- hy - ik y v i:X , (7) 
T n = (V x v n ) • z = ik x v n , y - ik y v n , x . (8) 

Note that in the reference frame in which k y = 0, Vi and T n are proportional to v> hy and 
v n>y , respectively. We combine Equations ([I])-© and after some algebraic manipulations we 
obtain the two following equations involving T; and T n only, namely 

Pi ~djT + a ' m ~dt pikz ° A 1 = a ' m ~dT' 

Pn^TT + air^n = 0> in T h (10) 

at 

where ca = B / JJTpi is the Alfven velocity. Note that the Alfven velocity is here defined 
using the density of the ionized fluid only. Equations (jHJ) and ( TlQjl are the governing equations 
for linear vorticity perturbations and, therefore, they are the governing equations of Alfven 
waves. 

For the subsequent analysis we define the ionization fraction, x, the ion-neutral collision 
frequency, u in , and the neutral- ion collision frequency, z/ ni , as follows 

Pn ^in *-^in /1 1 \ 

X = — , "in = , ^ni = • (11) 

Pi Pi Pn 

Since the collision frequencies are related by p x v m = p a u n i, we use v a i in all the follow- 
ing expressions for simplicity. Note that when p ; 7^ p n , v m 7^ u ni , meaning that the ion- 
neutral and neutral-ion c ollision frequencies are different (see a discussion on this issue in 



Zaqarashvili et al.l 1201 laf ) 



3. Normal Mode Analysis 

Here we perform a normal mode analysis. The temporal dependence of the perturbations 
is put proportional to exp (—icot), where u> is the angular frequency. From Equation (fTOj) we 
express T n in terms of F[ and insert the expression in Equation • We arrive at an equation 
involving T[ only, namely 

V(u,k K )Ti = 0, (12) 
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with 

V (u, k z ) = u 3 + i (1 + x) v ni uj 2 - k 2 z c\uo - iv ni k 2 z c\. 
For T; 7^ 0, the solutions to Equation ( TT2l) must satisfy V (uj, k z ) = 0, i.e. 

uj 3 + i (1 + x) VmU 2 - k 2 z c\uj - iv ni k 2 z c\ = 0. 



(13) 



(14) 



Equation (fl4|) is the dispersion relation of Alfven waves. Although with differen t notations 
Equat io n (fill) is equivalent to the dispersio n relati o ns previously f ound b y , e.g., 



Martin et al. 



( 1956 ): iKulsrud fc Pearcel Jl969h: iPudritzl (ll990h: 
jl998h : lKumar & Robert? fcp03h : IZaqarashvili et al. 
the absence of collisions, v n \ = and Equation f fT4|) becomes 



Piddington 



Jl997h: iKamava fc Nishi 



teOllbh : Mouschovias et all (120111 ). In 



w {u 2 - k 2 z c\) = 0. 



(15) 



From Equation (Tl5|) we get the classic dispersion relation of Alfven waves in an ideal plasma, 
namely u 2 = k 2 c\, and an additional mode with u> = 0. The general situation z/ n i ^ is 
investigated next. 



3.1. Standing waves 

We focus first on standing waves. Hence we assume a real wavenumber, k z , and solve the 
dispersion relation (Equation ( fl^|) ) to obtain the complex frequency, w = wr + icoi, with w R 
and cui the real and imaginary parts of cu, respectively. Since uj is complex the amplitude of 
perturbations is multiplied by the factor exp(co>it), with oj\ < 0. Therefore the perturbations 
are damped in time. 



Equation ([M]) is a cubic equation so it has three solutions. Unfortunately the exact 
analytic solution to Equation fTl4|) is too complicated to shed any light on the physics. 
However we can investigate the nature of the solutions using the concept of the polynomial 
discriminant. We perform the change of variable u = —is, so that Equation f lT4|) becomes 

s 3 + (1 + x) vms 2 + k 2 z c\s + u ni k 2 c 2 A = 0. (16) 

Equation ffl6|) is a cubic equation and all its c oefficients ar e real. From Equation ffT6]) we 



compute the discriminant, A, namely (see, e.g.. ICohenl 120001 ) 



A= - k 2 z c 2 A [A{l + xf^ 

- {x 2 + 20 X - 8) u 2 ni k 2 z c 2 A + 4ktc\] , (17) 

The discriminant, A, is defined so that (i) Equation ffl6|) has one real zero and two complex 
conjugate zeros when A < 0, (ii) Equation f fl6|) has a multiple zero and all the zeros are 
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real when A = 0, and (iii) Equation ffl~6]) has three distinct real zeros when A > 0. This 
classification is very relevant because the complex zeros of Equation (fl~6l) result in damped 
oscillatory solutions of Equation (fT4"|) whereas the real zeros of Equation ( fTS]) correspond to 
evanescent solutions of Equation (|14p . 



It is instructive to consider again the paradigmatic situation in which there are no 
collisions between the two fluids, so we set z/ ni = 0. The discriminant becomes A = —Ak®e\ < 
0, which means that Equation (1161) has one real zero and two complex conjugate zeros. 



Indeed, when z/ ni 



the zeros of Equation (TT6"j) are 

s = ±ik z CA, s 



which correspond to the following values of uj, 



±k z c A , 



to 



0. 



(18) 



(19) 



The two non-zero solutions correspond to the ideal Alfven frequency, as expected. 



We go back to the general case z/ n i ^ 0. To determine the location where the nature of 
the solutions changes we set A = and find the corresponding relation between the various 
parameters. For given u ni and % we fi n d two different values of k z , denoted by fc+ and k~, 
which satisfy A = 0, namely 



k± = — 



r 



20x - 8 ± \ ■ 



1/2 (X - 8) 3/2 
8(1 + X) 3 



-1/2 



(20) 



Since k z must be real, Equation (|20|) imposes a condition on the minimum value of \ which 
allows A = 0. This minimum value is x = 8 an d the corresponding critical k z is k^ = 
k~ = 3V3u ni /c A . When X > 8, Equation (EOJ) gives k+ < k z . For k z outside the interval 
(fc+, k~) we have A < so that there are two propagating Alfven waves and one evanescent 
solution. For k z £ (k^,k~) we have A > and so all three zeros of Equation ffT6l) are 
real, i.e., they correspond to purely imaginary solutions of Equation (1141) . There is no 
propagation of Alfven wa yes for k z £ (k^,k~). We call this interval the cut-off region. To 
the best of our knowledge, iKulsrud &: Pearcd (119691 ) were the first to report on the existence 
of a cut-off region of wavenumbers for Alfven waves in a partially ionized two-fluid plasma, 
when st udying the pro p agation of cosmic r ays. These cut-offs also appear in the works 



by, e.g., IPudritzl (11990f l: Kumar & Robertsl (120031 ) : Mouschovias et all (1201 If ) . The cut-off 
wavenumbers were ignored by IZaqarashvili et al.l (1201 lbl). who stated th a t there is always 
a solution of Equation ( IT4"1) with a non-zero real part. IZaqarashvili et al. ( 2011b ) probably 
reached this wrong conclusion because they never took x > 8 in their computations. Here 
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we clearly see that the three solutions of Equation 
and k z G k~). 



are purely imaginary when \ > 8 



Kamaya &: Nishil (119981 ) discussed the physical rea son for the existen ce of a range of 



cut-off wavenumbers in weakly ionized plasmas (see also iMouschoviasI 1 1 98 71 ) . When k z > k~ 
magnetic tension drives ions to oscillate almost freely, since the friction force is not strong 
enough to transfer significant inertia to neutrals. In this case, disturbances in the magnetic 
field affect only the ionized fluid as happens for classic Alfven waves in fully ionized plasmas. 
Conversely, when k z < kf the ion-neutral friction is efficient enough for neutrals to be 
nearly frozen into the magnetic field. After a perturbation, neutrals are dragged by ions 
almost instantly and both species oscillate together as a single fluid. The intermediate 
situation occurs when k z G {k^,k~). In this disturbance in the magnetic field 

decays due to friction before the ion-neutral coupling has had time to transfer the restoring 
properties of magnetic tension to the neutral fluid. In other words, neutral-ion collisions are 
efficient enough to dissipate perturbations in the magnetic field but, on the contrary, they 
are not efficient enough to transfer significant inertia to neutrals before the magnetic field 
perturbations have decayed. Hence, oscillations of the magnetic field are suppressed when 
k z G {k^,k~). Additional insight on the physical behavior of the perturbations near the 
cut-off region is given in Section 13.1.21 by analyzing the forces acting on the fluids. 

To avoid confusion we must inform the reader that the existence of cut-off wavenum- 
bers of Alfven waves discussed above is a purely two- fluid effect. These cut -off wavenum- 



bers are not the cu t -offs obtained in t he single-fluid appro xi mation (see, e.g.. iBalsaral Il996 



Forteza et al. 2007 ; Soler et al. 2009 ; Barcelo et al. 2011 ). Zaqarashvili et al. ( 20121 ) have 
shown that the cut-off wavenumbers found in the single-fluid approximation are a mathe- 
matical artifact, i.e., they are caused by the approximations made when proceeding from 
the multi-fluid equations to single-fluid equations and are not connected to any real physical 
process. On the contrary, the cut-off wavenumbers found in the two-fluid case are physi- 
cally and ma thematically real and are caused by the two-fluid interaction b etween ions and 



neutrals (see iPudritzl I199CH ; iKamaya &; Nishil 11998c iMouschovias et al.ll201ll ). 



3.1.1. Approximate analytic solutions 

We look for approximate analytic solutions to Equation (|T%j) corresponding to standing 
modes. We assume that k z is outside the cut-off interval (k^,k~), so that Equation ([Hj) 
has two complex solutions and one purely imaginary solution. This is the most interesting 
situation for the study of standing Alfven waves since no oscillatory modes exist when k z is 
within the cut-off interval. First we look for an approximate expression for the two oscillatory 
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solutions. To do so we write w = wr + and insert this expression in Equation ffT4"|) . We 
assume <C and neglect terms with oj\ and higher powers. Hence, it is crucial for 
the validity of this approximation that k z is not within or close to the cut-off region where 
% = 0. After some algebraic manipulations we derive approximate expressions for % and 
uj\. For simplicity we omit the intermediate steps and give the final expressions, namely 



2 



*?A + (1 + x? ^ 

X^ni , 2 2 



wi « ^ 2-^^4. (22) 

2[ky A + (i + x ) 2 u^ 

On the other hand, the remaining purely imaginary, i.e., evanescent, solution is u = ie, with 
the approximation to e given by 

m ^ci + (l + x)^' ( } 

When z/ ni = 0, we find Co>r, = ±A; 2 ca, wi = and e = 0, hence we recover the solutions in the 
uncoupled case (Equation (JT9"]) ). 

It is useful to investigate the behavior of the solutions in the various limits of v n \. First 
we consider the limit z/ ni <C k z cx, i.e., the case of low collision frequency, which means that 
the coupling between fluids is weak. Equations (j2"T|) and (123]) simplify to 

w R ~ ±k z c A , (24) 

wi « — , (25) 

e « -z/ ni . (26) 

In this limit % coincides with its value in the ideal, uncoupled case and oj\ is independent 
of k z . Hence, the the damping of Alfven waves does not depend on the wavenumber. 

On the other hand, when u n \ 3> k z CA, i.e., the case of strong coupling between fluids, 
we find 

- ±4=' ( 27 ) 

2^2 



X hie 



m « " x2 , (2* 



2(1 + X) 



e « -(l + xK- (29) 

Now the expression of cjr involves the factor i/l + x m the denominator, so that the larger 
the amount of neutrals, the lower cur compar ed to the value in the fully ionized case (see 



also lKumar fc Roberts! 120031 ; ISoler et al.ll2012l ). Now u\ is proportional to k%, meaning that 



the shorter the wavelength, the more efficient damping. 
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3.1.2. Comparison with numerical results 

Here we solve the full dispersion relation (Equation (I14p) numerically and compare the 
numerical solutions with the previous approximations (Equations f l2~Tj) - (r23]) ). First we set 
X = 2 and vary the ratio v^jk z CA between 10~ 2 and 10 2 . We compute o;r//c 2 ca and ui\jk z CA 
(see Figured]). The agreement between numerical and analytic results is very good. There is 
no cut-off region for this choice of parameters because we have taken x < 8- Regarding the 
real part of the frequency, we obtain that the oscillatory modes have wr/^ca ~ ±1 when 
Vin/k z CA <C 1. When the ratio v^ x jk z CA increases, un/k z CA decreases until de value wr/A^ca ~ 
±1/a/1 + x is reached. This behavior is consistent with the analytic Equation (12~TT) . The 
evanescent mode has wr = regardless the value of v ni jk z CA- The imaginary part of the 
frequency of the oscillatory modes tend to zero in both limits v ni jk z CA <C 1 and v^jk z CA 3> 1, 
while the damping is most efficient when v^jk z CA ~ 1- This result is also consistent with the 
analytic Equation f l22|) . although the approximation underestimates the actual damping rate 
when v n i/k z CA ~ 1- This discrepancy is a consequence of the weak damping approximation, 
which assumes |a;i| -C |u>r|. However when v n ^/k z CA ~ 1, the numerical results show that 
|co>i| and |cur| are of the same order, hence the damping is strong. The imaginary part of the 
frequency of the evanescent mode is very well approximated by Equation (|23|) . 

Next we increase ionization ratio to x = 20 and compute the same results as before 
(Figure [2]). Now there is a cut-off region because x > 8. The cut-off region is correctly 
described by Equation (T2Q|) . At the cut-off the oscillatory and evanescent modes interact 
and they become three purely imaginary solutions. Propagation is forbidden in this interval. 
As expected, the agreement between numerical and analytic results is not good near the 
cut-off region, but both results are in reasonably agreement far from the cut-off interval. 

To explore the physical behavior of the solutions near the cut-off region, we rewrite the 
momentum equations of ions (Equation ([1])) and neutrals (Equation (|2])) in the following 
forms, 

Pi 9 ^ = T-R, (30) 

Pn^ = R, (31) 

where T and R are the magnetic tension force and the friction force, respectively, given by 

T = -* Pi ^ Vi , (32) 

UJ 

R = Pn Vi- (33) 

UJ + 2Z/ ni 

In Equations fl30|) and fl3T|) we have not included magnetic pressure and gas pressure forces 
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Fig. 1. — Results for standing waves, (a) u^/k z CA and (b) uj\/k z c^ as functions of Vm/k z cx. 
We have used x = 1. Solid and dashed lines correspond to the numerical results of the 
oscillatory and evanescent modes, respectively, while the symbols correspond to the analytic 
expressions in the weak damping approximation (Equations ff2T|) (|23|) ) . 
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Fig. 2. — Same as Figure [U but with \ = 20. The shaded zone denotes the cut-off region 
according to Equation ( 1211 . 
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because they do not affect Alfven waves. Now we use the numerically obtained solutions for 
X = 20 (Figure [2]) to compute the moduli of T and R, namely ||T|| and ||R||, as functions 
of Vni/k z cx- Figure [3] displays the ratio ||T||/||R|| versus z/ ni //c^CA near the cut-off region. 
We have selected some locations in Figure [31 denoted by letters from a to e, to support the 
following discussion on the importance of the two forces. 

We start by analyzing the solutions on the left-hand side to the cut-off region. There 
are an oscillatory solution, a, and an evanescent solution, b. We find that ||T|| 3> ||R|| for 
the oscillatory solution a, so that there is a net restoring force for ions in Equation f)30p . 
Magnetic tension is the dominant force and drives ions to oscillate almost freely, whereas 
neutrals are only slightly perturbed by the weak friction force in Equation f l3Tj) . In the case 
of the evanescent solution b we obtain that ||T|| w ||R||. This means that there is no net 
force acting on ions. The evanescent solution b only produces perturbations in the neutral 
fluid. 

We turn to location c in Figure [3j i.e., within the cut-off interval. Here all the solutions 
are evanescent. The ratio | |T| |/| |R| | of the solution that was previously oscillatory decreases 
and becomes ||T||/||R|| < 1 before reaching the cut-off region. Now friction is the dominant 
force. Friction acts very efficiently in dissipating perturbations in the plasma before ions 
(and neutrals indirectly) have had time to feel the restoring force of magnetic tension. As a 
consequence, oscillatory modes are suppressed. 

Finally, we analyze the forces acting on the solutions on the right-hand side to the cut- 
off region. Again, there are an oscillatory solution, d, and an evanescent solution, e. As 
happened for the oscillatory solution a, we find that ||T|| > ||R|| for the oscillatory solution 
d, although in this case the friction force is not negligible. Magnetic tension provides now the 
necessary restoring force for the oscillations of ions, while the friction force is responsible for 
dragging neutrals when ions move. Hence, both species tend to oscillate together. Solution 
d represents a collective oscillation of the whole plasma. On the contrary, magnetic tension 
is negligible for the evanescent solution, e. This mode is governed by the friction force alone 
and simply causes the decay of perturbations. 

3.2. Propagating waves 

We move to the study of propagating waves. In the ideal fully ionized case the study of 
propagating waves is equivalent to that of standing waves. Here we shall see that ion-neutral 
collisions break this equivalence and propagating waves are worth being studied separately. 
For propagating waves we assume a real oj and solve the dispersion relation (Equation CHI) ) 
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to find the complex wavenumber, k z = fc 2i R + ik z j, with k Z} R and k Z) i the real and imaginary 
parts of k z , respectively. Then the amplitude of perturbations is multiplied by the factor 
exp(— k Z) iz), so that the perturbations are spatially damped. For real and positive u, k z ^ > 
corresponds to waves propagating towards the positive z-direction. Conversely, k z ,n < 
corresponds to waves propagating towards the negative ^-direction. In both situations, the 
sign of k z i is the same as that of k z ^. 



Equation ( 1T4"j) is a quadratic equation in k z . The solution to Equation (TT4"|) is 

U 2 U + i(l +X)Vni 



k , 



uj + iv ni 

= uj 2 jc\. For z/ ni / we write k z 



(34) 



When z/ ni = we recover the ideal result k 2 z = uj 2 jc\. For z/ ni / we write k z = k z ^ + ik z ^i 
and insert this expression in Equation fl34|) . Then, it is possible to obtain the exact expression 
of k 2 R , namely 



k 2 



ltU 2 L0 2 + (I + X)^ 



2c\ 



U 2 + V 2 



while the exact expression of k 2 j is 



k 2 



(^ 2 + (i + x)^i) 2 



u 2 u 2 + (i+ x y ni 



1/2 



(35) 



,R 



^ 2 + ^ 



(36) 



Contrary to the case of standing waves there is no cut-off region for propagating Alfven 
waves. It is always found that k z ^ ^ regardless the value of u. This apparent contradiction 
between the standing and propagating cases can be understood as follows. An oscillatory 
standing wave can be interpreted as the superposition of two propagating waves with the 
same frequency running in opposite directions. However such a representation does not work 
for a perturbation fixed in space and evanescent in time. The presence of static evanescent 
perturbations is a peculiar result only obtained when k z is real and uj is purely imaginary, 
a situation that cannot be described with propagating waves. For this reason evanescent 
solutions are absent from the present study of propagating waves. 



3.2.1. Approximate expressions and comparison with numerical results 



Going back to Equations fl35|) and fl36l) . we realize that it is possible to obtain simpler 
expressions of k 2 R and k 2 zl by taking advantage of the fact that the second term within 
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the square root of Equation (135]) is always much smaller than unity. We can approximate 
Equations f l35]) and (136]) as 

1 X 2zy ni^ 2 

fc 2 ~ 1 u " * 2 ^ 2 (38) 

As for standing waves we evaluate the analytic solutions (Equations (137]) and f ]38]) ) in the 
various limits of f ni . When ^ < w, Equations fl37]) and fl38]) simplify to 

« (40) 

Hence we find that k Zt & takes the same value as in the fully ionized case, and k z< i is indepen- 
dent of oj. Conversely when z/ ni ^> u we find 




4 4ci(l + X )Vf n 



« ±—y/l + x , (41) 
~ ± , XUj2 • (42) 

2(l + X )cA^in 

We find the presence of the factor i/l + x i n the approximate expression of k z ^ while now 
k Z) i is proportional to a; 2 . Hence, high-frequency waves are more efficiently damped than 
low-frequency waves. 

Now we solve the dispersion relation (Equation ffl~4]) ) numerically and compare the 
numerical solutions with the analytic approximations of k z ^ and k z j (Equations fl37]) and 
( 138]) ). These results and shown in Figure |4] for a particular set of parameters given in 
the caption of the Figure. As in the standing case, the agreement between numerical and 
analytic results is very good. We have replied these computations using various values of x 
(not shown here for simplicity) and the analytic results are always in accordance with the 
numerical ones. 
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3.2.2. Phase velocity and group velocity 

Here we compute the phase and group velocities of the propagating Alfven waves. The 
phase velocity, v p h, is defined as 

v P h = r4, (43) 

with ik denoting the unit vector in the direction of the wave vector. From the dispersion 
relation (Equation (fl4|) ) we get 



— = c A cos6M — - — . (44) 

k Y U + l(l + XK 

The phase velocity is a complex quantity. Its real part is related to the propagation speed of 
the wave while its imaginary part is related to the wave attenuation. Both real and imaginary 
parts of v ph have the same dependence on 9, meaning that the damping rate is independent 
of the direction of propagation. Ion-neutral collisions do not affect the dependence of v p h on 
9. From Equation ( I44p we can define the effective Alfven velocity, c"a, as 



U + Z(l + xKi 

The group velocity, v gr , is the propagation velocity of a wave packet and, therefore, of 
the wave energy. It is defined as 

_ duo A du A du A , . 

v - = VkW = wf* + w y ey + w z e - (46) 

Using the dispersion relation (Equation ( Fl4|) ) the group velocity can be more easily computed 

as 

'&D„ dV A dV \ 

e * + 7^r e v + i^r e z ' ( 4T ) 



that gives 



dV/dcu \dk x dk y y dk 

2k z c 2 A (ui + iu ni ) 
3u 2 - k 2 z c\ + i2(l + x)vmu 



0n fa + ^ni) 3/2 (U + i(l + X)^ni) V2 , s 

- ZCA 2^-2(i + x y m + i (i + x )u m u ez - m 

As it happens for the phase velocity, the group velocity is also complex. Since the damping 
of propagating Alfven waves depends on u, the components of the wave packet with higher u 
are more damped than the components with low u. Hence the concept of the group velocity 
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as the velocity at which the who le wave packet propagate s becomes obsolete when there 
is dissipation. However, following iMuschietti &: Duml ( 119931 ) it is still possible to define an 
effective group velocity as the velocity at which the center of the wave packet propagates. 
This effective group velocity is a time-dependent combination of the real and imaginary 
parts of v gr and depends on the particular form of the wave packet (see extensive details in 



Muschietti fc Duml 11993 ). 



In the limits of low and high collision frequency the damping is weak, so that the 
imaginary part of v gr can be neglected compared to the real part. When < u, we find that 
Vgr ~ cp^e z as in the ideal case, while when z/ ni 3> u, the group velocity is v gr ~ ca/ a/1 + X&z- 



Note again the presence of the factor yT + 



4. Initial-value problem 



In Section |3] we have followed a normal mode approach and have assumed that the 
temporal dependence of the perturbations is of the form exp (— iout). Here we go beyond the 
normal mode analysis and look for general time-dependent solutions to Equations and 
f fTOj) . To this end we follow two different approaches. First we solve the time dependent 
problem a nalytically by mea ns of the Laplace transform. Later we use the numerical code 
MolMHD ( iBona et al.l 120091 ) to numerically evolve the full set of basic Equations (pQ)— (jSJ- 
Finally, these two independent results are compared. In this section we r estrict ourselves to 
standi ng waves and assume a real and positive k z . We refer the reader to iRoberge fc Ciolek 
( 120071 ) for the problem of driven propagating waves in a two-fluid plasma. 



4.1. Analytic solution 



Here we look for general time-dependent solutions to the coupled Equations (J9j and ( TTUl) 
using the Laplace transform. We compute the Laplace transform in time of the vorticity 
perturbations T; and r n as 



T a (s) 



Tp(t) exp (—si) dt, 



(49) 



with s is the transformed variable and /3 denotes either T or 'n'. The Laplace transforms of 
the temporal derivatives of vorticity are given by 



C 



dt 



st p (s) - r 0l/3 , 



(50) 
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d 2 T{ 
dt 2 



/9> 



(51) 



where T 0j g and g denote the value of Tp and its temporal derivative at t — 0, respec- 
tively. From Equations ([9]) and ( TlOl) we deduce that the temporal derivatives of the vorticity 
perturbations at t — satisfy 



r 

1 r 



0,i 



V 



0,n 



-x^ni (r 0) ; - r 0)n ) , 
v ai (ro,i — r 0i n) . 



(52) 
(53) 



At the present stage, we leave r i and r n unspecified. 



We apply the Laplace transform to Equations ([9]) and (flOl) . From the transformed 
Equation (TlOl) we can express r n in terms of Y- x as 

1 



r n fs 



^niFi(s) + r , n 



(54) 



We insert this expression in the transformed Equation ([9]) and obtain the expression for Ti 

as 

fi( S ) = ^A, (55) 



where the functions B{s) and T>* 

B{s) = 

wis) = 



V*(s) 
(s) are defined as 



s [{s + v ni ) r ,i + x^r 0,nj ) 

s 3 + (1 + x) VmS 2 + k 2 z c\s + v ni k 2 z c 2 A , 



(56) 
(57) 



The function B(s) depends on the initial conditions and so it contains information about 
how the waves are excited. Conversely, the function T>*{s) is independent on the initial 
conditions. The function T>*(s) plays a very important role because it corresponds to the 
dispersion function. It tells us how the plasma behaves after the excitation has taken place. 
We note that the expression of T>*(s) (Equation (157|) ) coincides with the left-hand side of 
Equation (JTBI) . which was obtained from the normal mode dispersion relation (Equation (Till) ) 
after performing the change of variable u = —is. Hence the normal modes are consistently 
recovered from the present Laplace transform approach by setting D*(s) = 0. 

To compute the vorticity perturbation, T;, in the actual temporal domain we perform 
the inverse Laplace transform of Equation fl55l) . namely 



c- 1 




= c- 1 


\ B{a) 1 








_V*(s)_ 



(5? 



The inverse Laplace transform of Equati on (1581) can be evaluated by using the technique of 
partial fraction decomposition (see, e.g., Dykel 1999 ) if the roots of the equation V*(s) = 
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are known. The zeros of T>*(s) are precisely the normal modes studied in Section [31 so 
that we can take advantage of the approximate expressions found in Section [31 Once T^t) 
is known, the corresponding expression of r n (t) can be obtained using Equation (|54|) and 
applying the convolution theorem, namely 

F n (t) = r , n e-^ + u Qi [ r ; (r) exp [-u ni (t - r)] dr. (59) 

Jo 

From Equations and (15"9"|) we obtain T^t) and T n (t), namely 

Fi(t) = A 1 exp(et) 

+ [A 2 cos (tu R t) + A 3 sin (u R t)] exp (uit) , (60) 

T n (t) = Ci exp (et) + C 2 exp (-u ni t) 

+ [C 3 cos (ojRt) + C 4 sin (u;r£)] exp (u>it) , (61) 

with the approximate expressions of Co>r, Uj, and e given in Equations (l2Tj) - (l23]) and the 
coefficients A x -A 3 and C\-C± given in the Appendix. In the expression of wr (Equation (12T1) ) 
the + sign has to be taken. We must recall that the expressions of wr, wj, and e given in 
Equations f[2T]) fl23|) are computed in the weak damping approximation. 

In this analysis we have worked with the vorticity perturbations. However we can also 
write the solution to the initial-value problem using velocity perturbations. To do so we 
move, for convenience, to a reference frame in which k y = 0. Thus, the motions of Alfven 
waves are polarized in the y-direction and from Equations ([7]) and (|SJ) we get that the 
vorticity perturbations are proportional to the ^-component of velocities, namely v^ y and 
v D:y . Therefore the expressions for the temporal evolution of v- hy and the same as 

those given in Equations ( I6"0"j) and (IBTi) for T; and T n , respectively. We just need to replace 
To,i by and r 0i n by v n>0 , where f ij0 and v n>0 denote the values of the ion and neutral 
velocities at t — 0, respectively. 



4-1.1. Solution in the absence of magnetic field 



We start the study of some interesting cases by considering the situation in which there 
is no magnetic field . We set B = 0, and therefore ca = 0. This situation was studied by 
Vranjes et al.l (120081 ). Equations (J6"0j) and (J6l|) become 



+ 



r ,i + x r o,n 

i + x 
x (r ,i - r 0!ll ) 

i + x 



exp [- (1 + x) i/ni*] 



(62) 
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r n (t) 



r ,i + xTo, n 

i + x 
(Tpj — r 0n ) 

i + x 



exp [- (1 + x) ^mA ■ 



(63) 



These solutions agree with those found by IVranjes et al.l (120081 ). After an initial phase, i.e., 
when (1 + x) v ^ ^ 1) both r ; and T n relax to the same constant value, which corresponds 
to a weighted average of the initial conditions, namely 



Ti(t) « r n (t) 



To,i + xro,n _ Pir ,i + p n ro,n 



Pi + Pn 



(64) 



In the absence of magnetic field no oscillatory solutions are found. Vorticity perturbations 
remain constant after the relaxation phase. 



4-1.2. Low collision frequency 

Here we incorporate again the magnetic field but we consider the case of low collision 
frequency compared to the oscillation frequency, i.e., <C ur. Equations fl60l) and flBTl) 
simplify to 

Ti(t) « r 0>i cos (/c 2 c A t) exp (-^7^) > (65) 
r n (t) « r 0jI1 exp (-i/nit) . (66) 

When v n i <C cjr the vorticity perturbations of the ionized fluid and the neutral fluid are 
essentially decoupled from each other. After the excitation, vortical motions in the ionized 
fluid oscillate at the Alfven frequency, k z c\, and with a damping rate proportional to u ni . 
These are weakly damped Alfven waves. Conversely, vortical disturbances in the neutral 
fluid are evanescent in time. There is no driving force for vortical motions in the neutral 
fluid, and so there is no oscillatory behavior. 



4-1.3. High collision frequency 

Now we turn to the limit of high collision frequency compared to the oscillation fre- 
quency. This is the realistic case for many astrophysical applications and deserves special 
attention. We compute the limit of Equations (!60|) and (ISTj) when z/ ni 3> wr. We obtain 

r> u\ r °> i + * r °> n ( kzCA + 
Wt « ; COS — 7^=t 

1 + X \VTT^ 
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x exp I X 2 kzCA t 

' 1 2(1 + X ) 2 n 



+ X {T °:\ r °' n) exp [- (1 + X ) v m t] , (67) 
1 + X 

^ u\ r o,i + xr 0in / k z c A 
x e.pf ^_M£a, 

V 2(1 + X ) 2 "ni 

- r "'' ~ r "' n exp [- (1 + x) M . (68) 

By comparing Equations (IBTj) and ( 168]) with the equivalent solutions in the absence of 
magnetic field (Equations ( 162"]) and (163]) ) we find that in both cases there is a relaxation phase 
whose time scale is determined by the collision frequency. In the absence of magnetic field 
(Equations ( 162]) and ( 163]) ) the vorticity perturbations remain constant after the relaxation 
phase, while in the presence of magnetic field (Equations (IBTj) -( 168]) ) the perturbations of 
both ions and neutrals oscillate together as a single fluid at the modified Alfven frequency 
k z CA/ a/1 + X- The oscillations are exponentially damped in time. 

We define A(t) as the difference of and T n (t) computed from Equations (167]) and 
(168]) . namely 

A(t) = r ; (t) - r n (t) « (r 0ii - r , n ) exp (-—) , (69) 



Trel 



with r re i the relaxation time scale given by 



11 

Trcl = 7Z ~ v = TO 
(1 + X) "ni "in + "ni 

The function A(t) informs us about the relaxation phase. The function A(t) — > when 
t ^> r re i. Importantly, A(£) = when r 0j i = T Qn , i.e., there is no relaxation phase when the 
initial disturbance perturbs ions and neutrals in the same way. In such a case, both ions and 
neutrals oscillate together with amplitude equal to the initial condition. When r i ^ Y Q n 
the amplitude of the Alfvenic oscillations is the weighted average of the initial conditions 
(Equation (|64j)). The weights correspond to the densities of each fluid. 



4.2. Numerical solution 



Here we use the numerical code MolMHD to evol ye in time the bas ic Equations flTJ)- 
(]5]) in a fully numerical way. The numerical code (see iBona et all 12009k for details about 
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the scheme) uses the method of lines for the discretization of the variables, and the time 
and space variables are treated separately. For the temporal part, a 4th order Runge-Kutta 
method is used, whereas for the space discretization a finite difference scheme with a 4th 
order centered stencil is used. For a given spatial resolution, the time step is selected so as 
to satisfy the Courant condition. 

The MolMHD code evolves in time the plasma and magnetic field perturbations after 
an initial condition. In the present application the perturbations are assumed invariant in 
the x- and y-directions, so that perturbations only depend on the ^-direction and Alfven 
waves are strictly polarized in the ^/-direction. The only non-zero perturbations are the 
components of velocity and magnetic field in the ^-direction. The numerical integration of 
Equations ^T]) — (JSJ) is done in the interval z G [— 10L, 10L], where L is an arbitrary length 
scale. We use a uniform grid with 1001 grid points. Since we study standing waves the 
boundary conditions used at z = ±10L are that the velocity perturbations are fixed to zero. 
The boundary conditions for the other wave perturbations are that their spatial derivatives 
are set to zero. The initial condition at t = for the velocity of ions, v- hy , is 



which corresponds to the fundamental standing mode in our domain with dimensionless 
wavenumber k z L = n/20. In the fully ionized case, the dimensionless frequency of the 
fundamental mode is ujL/ca = 7r/20 ~ 0.157. Since we are dealing with linear perturbations 
we express the velocity perturbation in arbitrary units and set t>;.o = 1. 

Figure [5] shows the velocity perturbation of ions and neutrals at z = computed numer- 
ically with the MolMHD code for various values of the neutral-ion collision frequency. We 
compare the numerical results with the analytic ones (Equations (16T)|) and (RJTj) ). The results 
in the top row of Figure [5] (panels a, b, and c) are obtained when the velocity of neutrals 
at t = is the same as that of ions (Equation flTTl) ). while in the bottom row of Figure [5] 
(panels d, e, and f) neutrals are initially at rest. In the following paragraphs we discuss the 
results displayed in Figure Eland relate the time-dependent solutions with the normal modes 
investigated in Sectional 

(i) In Figure [5](a,d) we use z/ ni L/cA = 0.01 so that the collision frequency is an order of 
magnitude lower than the ideal Alfven frequency. Ions and neutrals behave almost indepen- 
dently. Regardless the initial condition for neutrals, ions oscillate at the Alfven frequency 
k z CA and are damped due to collisions. The oscillations of the ionized fluid act as a periodic 
forcing on neutrals, although the oscillations generated in the neutral fluid are of much lower 
amplitude compared to those in the ionized fluid. The dominant behavior in neutrals when 
the initial perturbation is nonzero is the exponential decay of the initial perturbation. This 




(71) 
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behavior is consistent with Equations f)65j) and fl66|) . In relation with the normal modes of 
Section [31 we see that when the collision frequency is low the oscillatory mode is mostly 
excited in the ionized fluid, while the evanescent mode dominates the neutral fluid behavior. 
Hence we can understand the physical reason for the existence of the evanescent mode as 
the way in which vorticity perturb ations in the neutral fluid decay in time due to collisions 
(see also IZaqarashvili et al.ll2011bl ). 

(ii) In Figure E(b,e) we increase the collision frequency to u ni L/cA = 0.1. The collision 
frequency and the ideal Alfven frequency are of the same order of magnitude. Now both 
ions and neutrals display strongly damped oscillations and it is not possible to relate the 
behavior of each fluid with one particular normal mode. Both oscillatory and evanescent 
normal modes are excited in both fluids and the observed behavior is the result of the joint 
effect of both oscillatory and evanescent modes. We notice that in Figure E^b,e) the analytic 
solution of the time- dependent problem does not exactly follow the full numerical result. The 
source of the discrepancy is in the approximate expressions of wr, u>i, and e (Equations (}2"T]) - 
( 123]) ) used in the analytic solution of the initial- value problem (Equations (1601) and (IBTj) ). 
These approximate expressions were obtained in the case of weak damping, which obviously 
does not apply to the situation of Figure E^b,e). However, it is remarkable that the analytic 
solution still captures the overall behavior correctly even beyond the weak damping limit. 
It is worth noting that instead of using the approximate values of cjr, ui, and e we could 
use their exact values obtained by numerically solving the dispersion relation. In such a 
case the agreement between analytic and numerical results is excellent (not shown here for 
simplicity) . 

(iii) Finally we use u^L/ca = 1 in Figure |5^c,f) so that the collision frequency is an order 
of magnitude higher than the ideal Alfven frequency. We find that ions and neutral behave 
as a single fluid and oscillate together at the modified Alfven frequency A; z ca/ a/1 + x- Again, 
both numerical and analytic results are in excellent agreement. The common amplitude of 
the ^-components of velocity of ions and neutrals, namely v y , is 

V y = . 72 

Pi + Pn 

The oscillations are weakly damped. When the initial velocity of neutrals is the same as 
that of ions (Figure [5](c)) both fluids oscillate together from t = 0. The evanescent mode 
is not excited when f i0 = f n ,o- Conversely, when neutrals are initially at rest (Figure E]T)), 
the single-fluid behavior begins after a short relaxation phase. Although very brief, the 
relaxation phase can be seen in Figure [5](f) near t = as sudden decrease of v\^ y and increase 
of v n ,y Figure [6] shows the same results displayed in Figure [5]T) but near t = in order to 
observe the relaxation phase in detail. The behavior of the velocity perturbations is governed 
by the evanescent normal mode during the relaxation phase and by the oscillatory normal 
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mode afterwards. This is consistent with Equations (1671) and fl68l) . In this case and since 
v a,o = 0) the amplitude of the joint oscillations after the relaxation phase (Equation f!72|) ) is 

Pi 1 , 7Q s 
v y = — ; %0 = — V i>0 . (73) 

Pi + Pn 1 + X 

Hence, the larger the ionization ratio \i the lower the amplitude. 



APPLICATION TO THE LOW SOLAR ATMOSPHERE 



In this Section we perform an application of the previous theoretical analysis to a partic- 
ular astrophysical plasma, namely the solar atmosphere. We specialize in solar atmospheric 
plasma because of two main reasons. The plasma in the coolest parts of the solar atmosphere, 
i.e., the photosphere and the low chromosphere, is weakly ionized. In addition, recent obser- 
vations have shown the ubiquitous presence of Alfvenic waves in th e solar atmosphere (e.g.., 



De Pontieu et al. 



2011 



2007 



To mczyk et al.l l2007t iMcIntosh et al.l 1201 ll ; lOkamoto fc De Pontieu 



De Pontieu et al.ll2012l ). It is therefore interesting to apply the theory developed in the 



waves 


due to ion- neutral collisions and its importance for plasma heating (see, e.g., 


Haerendel 


1992 




De Pontieu et al. 


2001; 


Khodachenko et al. 


2004; 


Leake et al. 


2005; 


Soler et al. 


2012|; 



Zaqarashvili et al.ll2013l . among others). In the present application we do not discuss damp- 
ing but investigate other effects on Alfven waves caused by partial ionization, namely the 
presence of cut-off wavelengths of standing waves and the modification of the effective Alfven 
velocity. For results about wave damping, the reader is referred to the papers cited above. 

We consider a simplified model for the solar atmospheric plasma. The variation of 
physical parameters with height from the photosp here to the base of the corona is taken from 
the quiet sun model C of IVernazza et al.l ( 19811 ). hereafter VALC model. To compute the 
friction coefficient, aa n , we ignore the influence of heavier species and consider only hydrogen. 
Hence, the expression of a- m is given in Equatio n (jSj), where we take a in ~ 10~ 20 m 2 according 
to t he estimation by IZaqarashvili et al.l (120131 ) for the case of direct elastic collisions (see 
also iBraginskiil Il965l ) . Figure [7] shows the inverse of the relaxation time of the ion-neutral 



coupling, t~2 (Equation (!72|) ). as a function of height in the low solar atmosphere. The 
very large values of r~j , i.e., very short r re i, point out that in the solar atmosphere ions 
and neutrals are very efficiently coupled and, for practical purposes, they can be considered 
as a single fluid. The single-fluid approach is usually followed in studies focused on wave 
damping (see the references in the previous paragraph). For comparison, we also plot in 
Fi gure [TJ the ion- neutral, y m , and neutral- ion, z/ ni , collision frequencies (see also similar plots 



m 



De Pontieu et al.ll200ll ). We find that the value of r re i is mainly determined by z/ in . Note 
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that in Figure [7] the solid and dotted lines are superimposed except at large heights. 



Concerning the magnetic field strength, we consider the model by lLeake fc Arberl (120061 ) 
of a vertical chromospheric magnetic flux tube expanding with height, whose form is given 
by 



B = B ph , (74) 

VPph/ 

where p = p\ + p n is the total density, B ph and p p h are the magnetic field strength and the 
total density, respectively, at the photospheric level, and /3 = 0.3 is an empirical exponent. 
We use pph = 2.74 x 10~ 4 kg m~ 3 from the VALC model and 2? p h = 1.5 kG. The dependence 
of p with height is determined by the VALC model. With this choice of parameters the 
magnetic field strength decreases with height so that B rj 100 G at 1,000 km and B 20 G 
at 2,000 km above the photosphere. 

The physical parameters in the simplified model of the low solar atmosphere used here 
depend on the vertical direction, z, while in the theoretical analysis of the previous Sections 
all the parameters are taken constant in z. However, it is possible to apply the expressions 
derived before for a homogeneous plasma to the present stratified case if the wavelength in 
the z-direction, \ z = 2n/k z , is much shorter than the length scale of the variation of the 
effective Alfven velocity. If such a restriction is fulfilled, we can perform a local analysis and 
use the physical parameters at a given height in the expressions derived for homogeneous 
plasma. We follow this approach. 

Another assumption made in this Section is that the amplitudes of the waves are small 
enough for the linear analysis to remain valid. This is a reasonable assumption for chromo- 
spheric waves. A parameter that can quantify nonlinearity of the waves is the ratio of the 
wave velocity amplitude to the local Alf ven velocity. The median of the velocity amplitudes 



of the chromospheric waves detected by lOkamoto fc De Pontieul ( 120111 ) is 7.4 km s 1 . This 



value is, at least, an order of magnitude lower than the expected value of the Alfven velocity 
in the chromosphere. 

We start by computing the cut-off region of wavenumbers, k z , of standing waves (Equa- 
tion (120]) ) as a function of height. Instead of k z , we show in Figure Ufa) the corresponding 
wavelength, X z . The shaded area in Figure E^a) indicates the range of wavelengths where 
oscillatory standing modes are not possible. At low heights above the photosphere the cut-off 
wavelengths are very short. The cut-off wavelengths increase with height and become of the 
order of a few kilometers at heights between 500 km and 1,500 km, approximately. Then 
the cut-off region disappears at a height of 1,600 km above the photosphere, approximately, 
because the ionization ratio, \i decreases with height as the plasma becomes more and more 
ionized and the threshold value x = 8 is reached at that height. 
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Next we turn to the computation of the effective Alfven velocity (Equation (|45j) ). 
Since this quantity depends on the wave frequency we use a frequency of 22 mHz, which 
corresponds to the dom i nant frequency of the chromospheric Alfvenic waves detected by 



Okamoto &: De Pontieul (120111 ). We display in Figure Efb) the real part of the effective 
Alfven velocity as a function of height. In Figure E(b) we also show the corresponding 
Alfven velocity computed using its classic definition that depends on the ion density only. 
We find that in the low chromosphere neutral-ion collisions are crucial to decrease the effec- 
tive Alfven velocity between one and two orders of magnitude with respect to the expected 
value taking into account the ion density only. 

Now we can check the assumption that X z is much shorter than the length scale of the 
variation of the effective Alfven velocity. On the one hand from Figure Efb) we get that 
the effective Alfven velocity varies about an order of magnitude in 2,000 km, approximately. 
On the other hand the largest cut-off wavelengths displayed in Figure Ufa) are of the order 
of a few tens of kilometers. Hence the use of a local analysis in the present application is 
justified. 



6. Implication for energy estimates 



In the solar context, the dissipation of Alfveni c waves may play an important role for 



De Pontieu et al. 2007 



Mcintosh et al. 



the h eating of the atmospheric plasma (see, e.g 
201ll ). The implications of part ial ionization f o r the calcul ations of en e rgy ca rried by Alfven 
waves were discussed by, e.g., IVranjes et al.l (120081 ) and iTsap et al.l ( 1201 II ). Both papers 



give expressions for the ene rgy flux of Alfven wa ves i n a partially ionize d plasmas. However 
the equations provided by Vranies et al.l (120081 ) and ITsap et al.l (120111 ) are different. The 



expression of I Vranies et a 



Eol) includes the factor (p-J p n ) , which is absent from the 



expression of ITsap et al.l ( 1201 ll ) . Hence the conclusions of these two papers regarding the 



i mpac t of partial ionization are in apparent contradiction. On the one hand, IVranjes et al. 



( 20081 ) explained that due to the factor (p;/p n ) the energy flux in weakly ionized plasmas 
beco mes orders o f magn itude smaller compared to the fully ionized case. On the other hand, 



since 



Tsap et al.l ( 120111 ) lacked that factor, they argued that the energy flux is independent 



of the ionization degree and obtained the same expression as for fully i onized plasmas . Here 
we s hall try to solve t his apparent contradiction between the results of IVranjes et al.l (120081 ) 
and ITsap et al.l (120111 ) . 



For application to the solar atmosphere we take the limit of high collision frequency. We 
choose a reference frame in which k y = 0. After the relaxation phase, i.e., for t ^> r re i, and 
neglecting the weak damping due to ion-neutral collisions, the ^-components of velocity of 



ions and neutrals oscillate together with frequency, k z cx/ a/1 + x? an d amplitude, v y , given in 
Equation (1721) . The energy flu x along magnetic field lines, S z , calculated using the Poynting 
vector (see, e.g., 



Walker 



20051) is 



S z — — — -£/ r &7/, 

2/i V1 



(75) 



where E x and 6^ denote the amplitudes of the x-component of the electric field and the 
y-component of the magnetic field perturbation. The electric field is 



E 



-Vi x B, 



(76) 



so that E x ~ v y B. On the other hand, the magnetic field perturbation is governed by 
Equation ([3]), from where we get b y ~ v y Byfl + x/°a- Hence the energy flux is 

B 



S z 



ca 2/x y 



2^ 



(77) 



with p = pi + p n the total plasma density. Consistently w e find the classical expression 
of the energy flux of Alfven waves (see, e.g., IWalkerl 120051 ) . with the only difference that 
the total plasma density replaces ion density and that the velocity amplitude after the 
relaxation phase, v y , is used. Since the relaxation phase is extremely short for realistic 
collision frequencies (Figure [7]), v y is the velocity amplitude that an observer would actually 
measure. In a fully ionized plasma, v y = v ii0 so that 

* = ir^ vl "- (T8) 



Up to here t here is no discrepancy between the analysis of IVranjes et al.l ( 120081 ) and 
Tsap et al.l ( ]201ll ). The discrepancy arises when the energy flux is written in terms of the 
initial velocity amplitudes of ions and neutrals. We insert Equation ( I72|) in Equation (177|) 
and rewrite the energy flux in terms of the initial velocity amplitudes, namely 

2^ V Pi + Pn J 

In weakly ionized plasmas p n ^> p x and the expression simplifies to 

B 



S* w tt^Vp ( — v i,o + w n ,o 

2^ VPn 



(80) 



In the a nalysis of IVranies et al.l (120081 ) the initial disturbance is localized in the ionized fluid 
only. So IVranjes et al.l (120081 ) took t> n;0 = 0. In this case Equation ( IHUl) becomes 
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Due to the presence of the factor (p-J p n ) the energy flux computed by lVranjes et al.l ( 20081 ) 
in a weakly ionized plasma is much lower than the corresponding value in a fully ionized 
plasma (Equation (|78p ). The physical reason is that significant porti on of the wave energ y 
is used to set into motion the neutral fluid, which is initially at rest in lVranjes et al.l ( 120081 ) . 



On the contrary, iTsap et al.l (120111 ) imposed v y = v- h o, which can be sat isfied only 
if fio = f n ,o according to Equation (!72|) . This means that ITsap et al.l (120111 ) implicitly 
assumed that both ions and neutrals are initially perturbed with the same velocity, although 
this condition is not explicitly stated in the paper. The energy flux is then 



B 



(82) 



The expression found bv lTsap et al.l (|2011) is t he same as Equation ( 178|) for a fully ionized 
plasma. This is so because in Tsap et al.l ( 2011 ) neutrals have initially t he same velocity as 
ions and no wave energy needs to be used in setting neutrals into motion. ITsap et al.l (120 111 ) 
claimed that the expression of IVranjes et al.l (120081 ) is incorrect. Here we clearly see that the 
discrepancy between Equations (IHTj) and (|82|) is caused by a different choice of the initial 
conditions and that both expressions are indeed correct. The source of the discrepancy was 
expressing the energy flux in terms of the initial velocity of ions. This is a velocity definitely 
not easy to measure observationally. Instead, the velocity that an observer can measure is 



the velocity amplitude after the relaxation phase, 
to estimate the energy flux. 



Then, Equation (ITTj) should be used 



Finally we recall that the present analysis, and so Equation ( 1771) . is valid in the case 
of Alfven waves propagating in a uniform medium so that the energy flux is uniform. As 
commented in Section [5j the application of Equation (ITT}) to the case of Alfvenic waves in the 
solar atmosphere must be done with caution si nce neither the densit y nor the wave velocity 
amplitude are uniform in the solar plasma (see iGoossens et al.ll2013l ). 



7. Conclusions 

In this paper we have studied theoretically Alfven waves in a plasma composed of an 
ion-electron fluid and a neutral fluid interacting by means of neutral-ion collisions. This 
configuration is relevant in many astrophysical and laboratory plasmas. To keep our in- 
vestigation as general as possible we have taken the neutral-ion collision frequency and the 
ionization degree as free parameters. 

First we have performed a normal mode analysis and have derived the dispersion relation 
of linear Alfven waves. The dispersion relation agrees with the equations previously found 
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( 120031 ) ; IZaqarashvili et al.l (l2011bf ); lMouschovias et al 



bv. e.g. . iKulsrud fc Pearce (11969 ): IPudritzl (119901): Martin et all (119971 ): iKumar fc Roberts 



(1201 if ). As in previous studies, we find 



that Alfven waves are damped by neutral-ion collisions. The damping is most efficient when 
the wave frequency and the collision frequency are of the same order of magnitude. This 
conclusion applies to both standing and propagating waves. The effective Alfven velocity of 
the plasma is modified by collisions, so that for high collision frequencies compared to the 
wave frequency the effective Alfven velocity depends on the t otal density of the plasma and 
not on the ion density only (see also lKumar fc Roberta 120031 ). 



In addition, an important result is that when \ > 8, i.e., when the plasma is weakly 
ionized, there is a range of wavenumbers (or, equivalently, of wavelengths) for which os- 
cillatory standing Alfven waves are not possible. These cut-off wavenumbers are physi- 



cal and are caused by a purely two- fluid effect (see IKulsrud fc Pearcd Il969l ; iPudritzl Il990 



Kamaya fc Nishilll998l : lMouschovias et al.ll201ll ). We investigated the physical reason for the 
existence of this cut-off region by analyzing the relative importance of the magnetic tension 
force and the neutral-ion friction force. We found that friction becomes the dominant force 
in the cut-off region, while tension is more important outside the cut-off interval. Hence, a 
disturbance in the magnetic field whose wavenumber is within the cut-off range decays due to 
collisions before the plasma is able to feel the restoring force of magnetic tension. As a conse- 
quence, oscillations of the magnetic field are effectively suppressed. The cut-off wavenumbers 
investigated here do not appear in the single-fluid approxim ation and are different f rom the 
unphysical cut-offs found in the single- fluid approximation (IZaqarashvili et al.ll2012f ). 



The solution of the initial-value problem is fully consistent with the normal modes 
analysis, and shows a growing strength of the interaction between ions and neutrals as the 
collision frequency increases. For high collision frequencies both ions and neutrals behave as 
a single fluid. Here a 'high collision frequency' means that the neutral-ion collision frequency 
is at least an order of magnitude higher than the wave frequency. This condition is fulfilled 
in many astrophysical plasmas which makes the single-fluid approximation appropriate in 
those situations for the computations of periods/wavelengths and damping times/damping 
lengths. However, as explained above the single-fluid limit does not fully capture the details 
of the interaction between ions and neutrals and, for example, the existence of a range of 
cut-off wavenumbers is a two-fluid result. 

As an example, we have considered Alfven waves in a plasma with physical condition 
akin to those in the low solar atmosphere. As a matter of fact, due to the large values 
of the collision frequency this plasma could be studied using the single-fluid approximation 
instead of the more general two- fluid theory. In that respect, the effective Alfven velocity 
is found to depend on the total density of the plasma. However, the presence of a certain 
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range of cut-off wavelengths that can constrain the existence of oscillatory standing modes 
in the chromosphere i s a result absen t from those previous studies that use the limit of 



strong coupling (e.g.. lHaerendell Il992l ; iDe Pontieu fc Haerendell Il998l ; iKhodachenko et al. 
2004 ; Leake et al. 20051 ). There are other astrophysical situations in which the ion- neutral 



coupling is not so strong as in the solar atmosphere and, therefore, a two-fluid treatment 
is needed. This may be the case of, e.g., protostellar discs (see a table with valu es of 
some physical parameters realistic of protostellar discs in iMalvshkin fc Zweibell 120111 ) and 



molecular clouds (see, e.g., lBalsaralll996t iMouschovias et al. 



20111 ). 



In the present work we have restricted ourselves to Alfven waves. Compressional mag- 
netoacoustic waves have not been investigated. It is expected that both types of MHD waves 
are simultaneously present in a plasma. Magnetoacoustic waves in a two-fluid plasma will 
be investigated in a forthcoming work. 
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22846 and from CAIB through the 'grups competitius' scheme and FEDER Funds. JT ac- 
knowledges support from MINECO through a Ramon y Cajal grant. RS thanks Ramon 
Oliver for reading the manuscript and for his constructive criticism. 
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A. Expressions of Coefficients 

The expressions of the coefficients of Equations fl60|) and (1ST]) are as follows, 

Ai = S € + V ^ + XV f^\ (Al) 
w| + (e - uji) 

a [^r + ^1 - e (u ni + 2m)} r ,i + x^mero.n , A2 ^ 
< + ( e - w l) 



This preprint was prepared with the A AS IATgX macros v5.2. 
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ujj (e + Vni + ui) - (Me- Ui) (z/ ni + m) r X"m (^R + ^I ~gg) r 

A 3 - — — 72 1 0,i H r 2 ^T i O,n, (Adj 



l|0/\2' \ / 

e + ^ni cu R + (e — wi) 
= t 1 + 7-7— ^71 r I J-V,,, (AO) 



K + Ki + Wi) ] (e 



°3 72 1 0,i - r 2 . ; ; ,21 r 9 . / 727 0,n, (AoJ 



w R + (e - z/ ni ) [o; R + (Vni + wi) ] + (e - wi) ] 

U R + ( e - ^ni) 2 WR [Wr + (fni + ^i) 2 ] [w R + (e - Wi) 2 ] 



with the expressions of wr, Wi, and e given in Equations (|2TT)-(l23l). In the expression of wr 
(Equation ( 12TT) ) the + sign has to be taken. 
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Fig. 3. — Ratio ||T||/||R|| versus v^ x jk z c^ for the solutions displayed in Figure [2] near the 
cut-off region (shaded zone). Solid and dashed lines correspond to oscillatory and evanescent 
solutions in time, respectively. 
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Fig. 4. — Results for propagating waves, (a) k z ^cx/uj and (b) k z ^c^/uj as functions of v^/uj. 
Solid lines correspond to the numerical results while symbols correspond to the analytic 
approximations (Equations (j3"7l) and ( |38i) ). We have used x — 2- 
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Fig. 5. — Results of the initial value problem. Velocity perturbation (in arbitrary units) 
of ions (solid line) and neutrals (dashed line) at z = as functions of time computed 
numerically with the MolMHD code for v^L/c^ = 0.01, 0.1, and 1. The symbols O and A 
are the corresponding analytic results (Equations f )60|) and fl6T]) ). In panels (a), (b), and (c) 
we have used the same initial condition for the velocity of both ions and neutrals, while in 
panels (d), (e), and (f) neutrals are initially at rest. We have used x = 2 in all cases. For 
reference, the dimensionless ideal Alfven frequency is uL/ca = vr/20 ~ 0.157. 
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Fig. 7. — Inverse of the relaxation time of the ion-neutral coupling (solid line) in the low solar 
atmosphere according to the VALC model. The ion-neutral (dotted line) and neutral-ion 
(dashed line) collision frequencies are shown for comparison. 
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Fig. 8. — (a) Cut-off region of wavelengths of standing Alfven waves as a function of height 
in the low solar atmosphere. The shaded area denotes the region where oscillatory standing 
modes are not possible, (b) Effective Alfven velocity (solid line) as a function of height for a 
wave frequency of 22 mHz. The Alfven velocity computed using the ion density only (dashed 
line) is shown for comparison. 



